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Abstract 

In this article we propose a new direct method of effective parameters restoration that is based on 
the wave propagation phenomenon. It retrieves the effective properties unambiguously, is applicable 
to thick metamaterial (MTM) slabs and is easy in implementation. It is validated on the case 
studies of fishnet, split cube in carcass, Jerusalem cross and ultrahigh refractive index MTMs. The 
constraints of the method are designated. 
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I. INTRODUCTION 



MTMs research is a new emerging field of pfiotonics due to extravagant properties tfiat 
can be obtained after relevant engineering of their structure. It provides great opportunities 
to control the light propagation. The potential applications of MTMs range from invisibility 
cloaksi to superresolution lenses^ and nano couplers^, but are not limited to them. 

After some successful proof-of-principle experiments with the negative index behaviour, 
MTMs are moving now towards new frontiers, for example, having permittivity e or per- 
meability /i either close to zero or very large, exhibiting isotropy or extreme anisotropy, 
dumping or enhancing spatial dispersion. In the cases when effective parameters (refractive 
index n, impedance Z, permittivity e and permeability yu) can be introduced in the correct 
way, reliable retrieval procedure is on demand. Assigning some certain values for the param- 
eters is the mean by which the properties of the MTM can be characterized in the universal 
mode. 

The existing retrieval methods can be divided into following groups: 

1) Reflection/transmission based. These methods are based on the reversal of the ampli- 
tude and phase reflection r and transmission t (or S-parameters) from the slab of homogenous 
material. Initially proposed for normal incidence on the isotropic metamaterial^*^ the pro- 
cedure was modified for inclined incidence-, bianisotropic^'*^ and chiral metamaterials^. The 
procedure introduced by Smith et al.- is easy to implement and it is widely used, so let we 
call it the standard method (SM). 

SM suffers from the well-known "branch" problem. As we restore refractive index Re{n) 
using the arccos function we have ambiguity of 27rm/(/c(i), where m is an integer number, k 
- wavenumber and d - slab thickness. This method works well only for thin slabs. Unfor- 
tunately, a thin slab consisting from just a few MTM layers exhibits the properties which 
can be very different from the bulk MTM—. Applying standard method to thick MTM slab 
is usually problematic. First, different branches of Re{n) become very close to each other. 
Second, MTMs are usually absorptive so transmission through the thick slab is very low 
that leads through a numerical noise to large errors. 

2) Field averagingiiii^. In this case the fields {E, D, H, B) in the unit cell are averaged on 
some surfaces and lines and the effective properties are found according to material relations 
D = eE and B = fiH. Determining appropriate surfaces and contours for averaging is a 
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non-trivial problem especially when optimizing 3D MTM structures. 

3) Recently published quasimode theoryi^. This method is based on the maximization of 
optical density of states for metamaterial while changing e and /i of the surrounding medium. 
The method is computationally demanding as it requires 4-parameters optimization for each 
frequency. 

4) Based on the wave phenomenon. In this case the effective parameters are derived 
from observation of the wave propagation phenomenon. For example, Popa and Cummer- 
proposed to match the simulated fields inside and outside of the thin metamaterial slab with 
the fields theoretically predicted by formulae. Fitting parameters in this case are refractive 
index and impedance. Due to multiple reflections the standing wave forms inside the slab 
so the fitting procedure is complicated. 

The method we propose in this article is also based on the wave propagation phenomenon 
and we call it wave propagation retrieval method (WPRM). In some sense it exploits and 
elaborates ideas proposed in Ref.= but WPRM has significant difference. We propose to 
restore the refractive index n from the wave propagation phenomenon through very thick 
metamaterial slab. In this case we deal with only one wave propagating in the metamaterial 
so the retrieval formulae become easy to use. Such retrieval is unambiguous, and the effective 
parameters derived from the thick layer are bulk automatically. Impedance z is restored from 
the input interface reflection only that leads to uniqueness of the retrieval and absence of 
effects caused by flniteness of the structure. 

The article is organised as follow. In the part [III we give mathematical foundations of 
WPRM and systematically deflne the limits of its application. Case studies demonstrate 
the power of WPRM in the part lllli Conclusions sum up the article. 

II. WAVE PROPAGATION RETRIEVAL METHOD 
A. Methodology 

The simulation geometry assumes the MTM slab infinite in x— and directions and 
semi-infinite in direction. The incidence of the wave occurs from vacuum to MTM. Inside 
MTM there is only one plane wave propagating along 2;— axis (Fig. [11(a)). The fields have 
complicated spatial distribution inside the unit cell, so in order to apply an effective medium 
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FIG. 1: Wave propagation simulation configuration (a); average amplitude and phase dependence 
on the unit cell number (theoretical expectation) (b). 

description we should abstract out from composite internal structure and regard it as a bulk 
homogeneous material. Thus we average the components of complex field vectors over the 
unit cell. Our choice is the electric field, but the procedure can be easily reformulated for 
the magnetic field. One more assumption we utilize here for the sake of simplicity is the 
linear polarization of the fields, which lead to the scalar procedure of restoration. 

As we are working in the effective medium approach we describe the light wave of fre- 
quency u; as a modulated plane wave: 

E{x, y, z) = Eou{x, y, z)exp{—ik'z — k"z), (1) 

where u{x,y,z) is a modulation function periodic in z— direction with the same period a 
as MTM, k' = n'u/c and k" = n"uj/c are real and imaginary part of the wavenumber 
connected with the real n' and imaginary n" parts of the refractive index. Modulation 
function is normalized over the unit cell volume V: J u{x,y, z)dxdydz/V = 1. It changes 
quickly over the unit cell comparing with the exponential ansatz in ([1]). Assuming the slow 
variation of the phase, which is valid under the generally accepted homogeneity conditions 
Xeff > 4:a, the averaged over the A^— th unit cell electric field is equal to 

{E)j^ = J EQu{x,y, z)exp{—ik' z—k" z)dxdydz ~ EQexp{—k" ZN)exp{—ik' zn) = Aexp^—icp), 

(2) 

where zn = {N — l/2)a is the coordinate of the A^— th unit cell centre, A = \{E)m\ and 
= arg((£')7v) are the averaged field amplitude and phase, and the integration is performed 
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over unit cell's volume. 

Both parts of the refractive index are straightforwardly derived via integral changes in 
phase and amplitude logarithm A In A over the fixed number of unit cells AN: 

A0 c 



n 



AN au 



(3) 



"=-^^- 

Semi-infinite in 2;— direction configuration means that there is only one interface between 
vacuum and MTM. In this case impedance Z can be restored in the unique way from the 
reflection r: 

Z = — . (5) 

It was mentioned by Landau^ that knowledge of the impedance for absorbing medium 
allows determining the fields outside it. In our case we solve the reverse problem determining 
impedance from the external fields. Knowing the refractive index and impedance one can 
easily obtain permittivity and permeability e = n/ Z, ^ = nZ . 

We performed the simulations of the wave propagation inside the thick MTM slab (20- 
50 MTM layers) with CST Microwave Studio 2009 (CST Computer Simulation Technology 
AG, Darmstadt, Germany) with a time-domain solver. For broadband excitation short 
modulated Gaussian pulse was chosen. As a reference we used the effective parameters 
retrieved with a SM from the reflection/transmission through 3 MTM layers simulated with 
a frequency-domain solver in CST. We found out that results retrieved with WPRM are the 
same as retrieved with SM for the pulse relative bandwidth up to 100% (see part IIII DP . 
Averaging and retrieval according to the formulae were made with a homemade Matlab 
program. An example of WPRM realization for a split-cube in carcass (SCiC) structure 
(design description is below in the text) is shown on the Figure [2j 



B. WPRM advantages and limitations 

WPRM has several clear advantages. First, it does not experience the "branch" problem. 
Second, WPRM does not suffer from low transmission and numerical noise as it uses the 
amplitude and phase distribution inside the structure - not only on the first and last slab 
interfaces. The method works well for thick slabs. Third, restoring the effective parameters 
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FIG. 2: An example of (a) amplitude and (b) phase of the averaged field distribution over 20 unit 
cells of the Split-Cube-in-Carcass design for the negative refractive index regime. Colour scales 
used in both cases for demonstration of the averaged amplitudes and phases are not equivalent. 



with SM one should always use the broadband excitation source. There is no difference for 
WPRM whether to use broadband or narrowband source, so one can choose what excitation 
is more appropriate for a given problem. 

WPRM also has some limitations. We cannot simulate semi-infinite MTM in practice so 
the reflection from the back interface may happen and lead to the standing wave formation. 
Different amplitude In A distributions can occur except the most desirable linear one (Fig. 
[3t^a)) when formula (jlj) works well. In the case of low wave attenuation and reflection from 
the back interface we can observe a standing wave pattern (Fig. [3](b)). To retrieve n and Z 
correctly we should make the MTM slab thicker and take into account only the first linear 
part of the graph. In the case of large wave attenuation (Fig. [3](c)) the amplitude value 
drops down within a short distance to the level of numerical noise. Noise data should be 
excluded for restoration. 

To find out the limitations of WPRM we simulated "the worst homogenous 
slab with pre-defined refractive index n and impedance Z = 1 with a perfect mirror (PEC 
boundary condition) on the other side of the slab (Fig. |4](a)). Numerical simulations were 
performed in CST until the total electromagnetic energy stored in the system dropped down 
to -80dB with respect to the maximum level. 

It is obvious that not the value of the refractive index itself but the terms in the ex- 
ponentials in ([2]) play the role in phase and amplitude changes. So we calculated relative 
error of n restoration for different n'kod and n^k^d (Fig. Ill^b)-(d)), where = u/c is the 
wavenumber in vacuum. In the case of fast amplitude decay (Fig. [3](c)) we should consider 
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FIG. 3: Possible cases of simulated average amplitude distribution inside the MTM slab: (a) normal 
case; (b) standing wave formation due to the back interface reflection and slow wave attenuation; 
(c) large wave attenuation and numerical noise. 

d not the full thickness of the slab but only until the In A reaches the noise level. Thou 
different regions correspond to the small restoration error in n' and n" (Fig. H] (b) and (c)), 
there is an optimal region where the total relative error is less than 1% (Fig. |4](d)), which 
can be approximated with a rectangle. We checked WPRM for positive index materials 
but the results are applicable for the negative index materials as well. So we might state 
that WPRM is applicable with a relative error less than 1% for —45 < n'kod < 45 and 
5 < n"kod < 56. In the usual MTM cases these limits are wider as we do not have perfect 
reflector at the second interface. 

The range of applicability of WPRM can be extended to the case of low-loss structures 
if we excite the MTM slab with short Gaussian pulse and stop the simulation when the 
centre of the pulse reaches the second interface of the slab. In this case we have a reflected 
wave only near the second interface and can skip these amplitude and phase data. As 
a shortcoming such simulation is possible only for metamaterials with low group velocity 
dispersion and requires additional estimation of the time the pulse travels through the slab. 
The case studies of applicability of WPRM to low-loss MTMs follow below in the sections 
ImCl and UTTDl 

III. CASE STUDIES 

To check how good WPRM works we choose MTM designs for different electromagnetic 
range: optical, terahertz and microwave. Standard method^ was used as a reference. 
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FIG. 4: Configuration of "the worst case" simulation (a). Relative error of (b) n', (c) n" restoration 
and (d) total relative error. Black area correspond to the relative error less than 1%. 

A. Fishnet 

Fishnet structure^^ is a classical example of the optical negative index MTM with high 
transmittivity. For WPRM test we took a specific design parameters from^°. It was shown 
iii^ the dominance of a single Bloch mode and fast convergence of the effective parameters 
with a slab thickness. Effective parameters n and Z restored with WPRM and SM are in a 
perfect agreement (Fig. [5]). 

B. Split cube in carcass 

The split cube in carcass structure (SCiC) is similar to the split cube in cage desigr>i^. 
Its remarkable property is bulkiness: the effective parameters at negative index region do 
not depend on the number of layers starting from a single layer—. The structure possesses 
cubic symmetry that insures its polarization insensitivity. As a drawback, SCiC has quite 
low transmittivity and high absorption. 

The SCiC unit cell consists of two silver parts embedded in silica (refractive index n = 
1.5). Silver is regarded as Drude metal with plasma frequency Up = 1.37 x 10^^ rad/s and 
collision frequency = 8.50 x lO^^s"^— . The outer part carcass (Fig. [H] (a)) is a sort of 3D 
wire medium and provides negative e'. The inner part split cube (Fig. E] (b)) is the hollow 



8 



(a) p 

1.5 - 




160 



180 200 
Frequency (THz) 



220 



(b) 




N 

O 

u 
c 

(S 



0.5 



Re(Z) WPRM 

Re(Z)SM 

- - - lm(Z) WPRM 
~ " - lm(Z) SM 



L 







-t \ — 




160 180 200 220 
Frequency (THz) 



FIG. 5: Refractive index (a) and impedance (b) of the fishnet structure restored with WPRM 
(red hne) and SM (green hne). Real and imaginary parts are marked with sohd and dashed hnes 
respectively. This colour and symbol convention will be respected in all subsequent figures. 

cube with the slits in the middle of the facets. It is the logical 3D extension of the split ring 
resonator concept and gives negative fi'. The details about structure's sizes are labeled on 
the Figure [61 

Both WPRM and SM give very similar results (Fig. [7]). For the frequencies above 190THz 
there are differences in the retrieved impedance (Fig. [7] (b)). They come out from the fact 
of the dipole resonances of split cubes and failure of the effective homogeneity approach^. 

C. Jerusalem crosses 

Recently Wang et al. reported a Jerusalem crosses 3D MTM design for microwave region 
claiming isotropic broadband left-handednesa^. As the reported losses in this MTM are 
relatively low it was a good opportunity to check WPRM. All the design parameters were 
taken the same as ir>^. 

To our astonishment WPRM showed no negative refractive index in the specified region 
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FIG. 6: Split Cube in Carcass unit cell design: (a) Carcass, (b) Split Cube. 
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FIG. 7: Split Cube in Carcass effective parameters: (a) refractive index and (b) impedance. 

(around 6GHz) (Fig. [8]^a)). WPRM results perfectly coincide with the SM results if one 
chooses the correct branch. These results are feasible as real and imaginary parts of n must 
obey the Kramers-Kronig relations^^. Small imaginary part corresponds to smooth changes 
in the real part in contrary tc^^. This is another confirmation that WPRM allows to avoid 
errors and misunderstanding in the refractive index restoration and interpretation. 

A kink in the impedance spectrum restored with SM (Fig. [8]^b)) appears due to the 
finiteness of the system. Additional studies showed that it changes its position with the 
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FIG. 8: Jerusalem crosses MTM effective parameters: (a) refractive index and (b) impedance. 

number of layers always being exactly at the frequency of Fabry-Perot resonance. There is 
no such problem for WPRM as multiple reflections from the slab interfaces are successfully 
avoided. 



D. 3D metamaterial with ultrahigh refractive index 

A 3D MTM with ultrahigh refractive index (URIMTM) over a broad bandwidth has been 
recently proposed by Shin et al.^''. Their idea is to use unit cell consisting of metallic plates 
connected with wires to obtain the refractive index up to 7. In^S there were not specified the 
exact geometrical parameters of the unit cell so we constructed a similar structure (Fig. [9]) 
from silver (metal parameters the same as in part IIIIBj) embedded in silicon {e = 12). 



URIMTM shows high refractive index up to 7 in a broad range of frequencies (Fig.[TO](a)). 
WPRM and SM results are very close to each other. A kink in the impedance around 30 
THz (Fig. [TUlfb)) is due to finiteness of the structure as it is discusses in the part IIII CI 

It is worth to mention that lg{n"kod) = —0.75 at frequency 30THz, which is out of the 
low error limits defined in Figure Hl^d), but relative difference of WPRM with SM is only 
7%. This fact confirms applicability of WPRM to low-loss MTM as well as to low-loss. 
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FIG. 9: URIMTM design geometrical parameters. Only half of the unit cell is shown. Metallic 
structure is embedded in silicon. 
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FIG. 10: URIMTM effective parameters: (a) refractive index and (b) impedance. 
IV. CONCLUSIONS 



In this work we proposed a new effective parameters retrieval method based on the wave 
propagation phenomenon and validate its applicability to lossy and low-loss metamaterials 
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on the case studies. The main advantage of WPRM is that it retrieves the effective properties 
unambiguously and is apphcable to a thick MTM slab. WPRM does not depend on the 
working frequency so it can be used in microwaves, terahertz and optical regions - while an 
effective medium approximation is valid. WPRM is easy to implement using commercially 
available electromagnetic simulation software. 

In principle, WPRM can be extended to anisotropic media or media with a non-local 
response as soon as required observable features of the wave can be detected numerically. 
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